####################################################################
## author:    Robert A. Huber
## contact:   robert.huber@ir.gess.ethz.ch
## file name: pc_uk_boot.R
## Context:   Populism and Climate Sceptism, individuals from BES
## started:   2016-10-12
## Summary:   bootstraps factors as RBC
######################################################################

df_bes$populism1_num <- ifelse(df_bes$populism1 == "Don't know", NA,
                                as.numeric(df_bes$populism1))
df_bes$populism2_num <- ifelse(df_bes$populism2 == "Don't know", NA,
                                as.numeric(df_bes$populism2))
df_bes$populism4_num <- ifelse(df_bes$populism4 == "Don't know", NA,
                               as.numeric(df_bes$populism4))
df_bes$populism5_num <- ifelse(df_bes$populism5 == "Don't know", NA,
                               as.numeric(df_bes$populism5))
df_bes$populism6_num <- ifelse(df_bes$populism6 == "Don't know", NA,
                               as.numeric(df_bes$populism6))

df_bes$efficacy1_num <- ifelse(df_bes$efficacyPolCare == "Don't know", NA,
                                             as.numeric(df_bes$efficacyPolCare))

df_bes$efficacy2_num <- ifelse(df_bes$efficacyNoMatter == "Don't know", NA,
                                              as.numeric(df_bes$efficacyNoMatter))

df_bes$efficacy3_num <- ifelse(df_bes$efficacyNotUnderstand == "Don't know", NA,
                                                   as.numeric(df_bes$efficacyNotUnderstand))

df_bes$efficacy4_num <- ifelse(df_bes$efficacyUnderstand == "Don't know", NA,
                                                -as.numeric(df_bes$efficacyUnderstand))

df_bes$efficacy5_num <- ifelse(df_bes$efficacyTooMuchEffort == "Don't know", NA,
                                                   as.numeric(df_bes$efficacyTooMuchEffort))

df_bes$auth1_num <- ifelse(df_bes$auth1 == "Don't know", NA,
                               as.numeric(df_bes$auth1)-1)

df_bes$auth2_num <- ifelse(df_bes$auth2 == "Don't know", NA,
                           as.numeric(df_bes$auth2)-1)

df_bes$auth3_num <- ifelse(df_bes$auth3 == "Don't know", NA,
                           as.numeric(df_bes$auth3)-1)

df_bes$auth4_num <- ifelse(df_bes$auth4 == "Don't know", NA,
                           as.numeric(df_bes$auth4)-1)

fut_varspop <- paste0("populism", c(1:2, 4:6), "_num")
fut_varseff <- paste0("efficacy", c(1:5), "_num")
fut_varsaut <- paste0("auth", c(1:4), "_num")

nboot <- 1000

cUK <- " + auth + gender + inc + age  + edu_high + riskTaking + polAttention + satDem + efficacy + genTrust + econPersonalRetro +  econGenRetro  + partyId_ukip + partyId_con + partyId_lab + partyId_lib + partyId_gre + partyId_oth"  # object with "fixed control variables"
outCC1 <- matrix(NA, nrow = nboot, ncol = 3)
outEG1 <- matrix(NA, nrow = nboot, ncol = 1)
outEP1 <- matrix(NA, nrow = nboot, ncol = 3)
for (i in 1:nboot) {
  
  bootdf <- df_bes[sample(nrow(df_bes), replace = TRUE),]
  bootdf$populism <- principal(bootdf[, fut_varspop], nfactors=1, rotate="oblimin")$scores[,1]
  bootdf$efficacy <- principal(bootdf[, fut_varseff], nfactors=1, rotate="oblimin")$scores[,1]
  bootdf$auth <- principal(bootdf[, fut_varsaut], nfactors=1, rotate="oblimin")$scores[,1]
  
  boot_outCC1 <- multinom(formula(paste0("climateChange ~ populism + lr" , cUK)), data = bootdf, weights = bootdf$wt_full_w7, maxit = 1000)
  boot_outEG1 <- lm(formula(paste0("enviroGrowth ~ populism + lr", cUK)), data = bootdf, weights = bootdf$wt_full_w7)
  boot_outEP1 <- multinom(formula(paste0("enviroProtection_RBC ~ populism + lr", cUK)), data = bootdf, weights = bootdf$wt_full_w7, maxit = 1000)
  
  outCC1[i,] <- coef(boot_outCC1)[,"populism"] # n X k matrix of coefficients
  outEG1[i,] <- coef(boot_outEG1)["populism"] # n X k matrix of coefficients
  outEP1[i,] <- coef(boot_outEP1)[,"populism"] # n X k matrix of coefficients
  
  #pred_out_min[i] <- predict(boot_out, newdata = blah) # n X 1 vector of predictions
  #pred_out_max[i] <- predict(boot_out, newdata = blah2) # n X 1 vector of predictions
  
  #pred_out$y11[i] <- predict()
  #pred_out$y10[i] <- predict()
}

quantile(outCC1[,1], c(0.05,.5, 0.95))
quantile(outCC1[,2], c(0.05,.5, 0.95))
quantile(outCC1[,3], c(0.05,.5, 0.95))

quantile(outEP1[,1], c(0.05,.5, 0.95))
quantile(outEP1[,2], c(0.05,.5, 0.95))
quantile(outEP1[,3], c(0.05,.5, 0.95))

quantile(outEG1[,1], c(0.05,.5, 0.95))

write.csv(data.frame(outCC1), file = paste0("./boot data/outCC1_", Sys.Date(), ".csv", sep=""))
write.csv(data.frame(outEP1), file = paste0("./boot data/outEP1_", Sys.Date(), ".csv", sep=""))
write.csv(data.frame(outEG1), file = paste0("./boot data/outEG1_", Sys.Date(), ".csv", sep=""))

varList <- c("populism", "lr", "populism:lr")

outCC2pop <- matrix(NA, nrow = nboot, ncol = 3)
outEG2pop <- matrix(NA, nrow = nboot, ncol = 1)
outEP2pop <- matrix(NA, nrow = nboot, ncol = 3)
outCC2poplr <- matrix(NA, nrow = nboot, ncol = 3)
outEG2poplr <- matrix(NA, nrow = nboot, ncol = 1)
outEP2poplr <- matrix(NA, nrow = nboot, ncol = 3)

for (i in 1:nboot) {
  
  bootdf <- df_bes[sample(nrow(df_bes), replace = TRUE),]
  bootdf$populism <- principal(bootdf[, fut_varspop], nfactors=1, rotate="oblimin")$scores[,1]
  bootdf$efficacy <- principal(bootdf[, fut_varseff], nfactors=1, rotate="oblimin")$scores[,1]
  bootdf$auth <- principal(bootdf[, fut_varsaut], nfactors=1, rotate="oblimin")$scores[,1]
  
  boot_outCC2 <- multinom(formula(paste0("climateChange ~ populism * lr" , cUK)), data = bootdf, weights = bootdf$wt_full_w7, maxit = 1000)
  boot_outEG2 <- lm(formula(paste0("enviroGrowth ~ populism * lr", cUK)), data = bootdf, weights = bootdf$wt_full_w7)
  boot_outEP2 <- multinom(formula(paste0("enviroProtection_RBC ~ populism * lr", cUK)), data = bootdf, weights = bootdf$wt_full_w7, maxit = 1000)
  
  outCC2pop[i,] <- coef(boot_outCC2)[,"populism"] # n X k matrix of coefficients
  outEG2pop[i,] <- coef(boot_outEG2)["populism"] # n X k matrix of coefficients
  outEP2pop[i,] <- coef(boot_outEP2)[,"populism"] # n X k matrix of coefficients
  
  outCC2poplr[i,] <- coef(boot_outCC2)[,"populism:lr"] # n X k matrix of coefficients
  outEG2poplr[i,] <- coef(boot_outEG2)["populism:lr"] # n X k matrix of coefficients
  outEP2poplr[i,] <- coef(boot_outEP2)[,"populism:lr"] # n X k matrix of coefficients
  
}


write.csv(data.frame(outCC2pop), file = paste0("./boot data/outCC2pop_", Sys.Date(), ".csv", sep=""))
write.csv(data.frame(outEP2pop), file = paste0("./boot data/outEP2pop_", Sys.Date(), ".csv", sep=""))
write.csv(data.frame(outEG2pop), file = paste0("./boot data/outEG2pop_", Sys.Date(), ".csv", sep=""))
write.csv(data.frame(outCC2poplr), file = paste0("./boot data/outCC2poplr_", Sys.Date(), ".csv", sep=""))
write.csv(data.frame(outEP2poplr), file = paste0("./boot data/outEP2poplr_", Sys.Date(), ".csv", sep=""))
write.csv(data.frame(outEG2poplr), file = paste0("./boot data/outEG2poplr_", Sys.Date(), ".csv", sep=""))

#### IF BOOT RAN ALREADY ####

nboot <- 1000

outCC1 <- read.csv("./boot data/outCC1_2018-05-15.csv")
outCC1 <- outCC1[,2:4]
outEG1 <- read.csv("./boot data/outEG1_2018-05-15.csv")
outEG1 <- outEG1[,2]
outEP1 <- read.csv("./boot data/outEP1_2018-05-15.csv", colClasses = c(NA, "double", "double", "double"))
outEP1 <- outEP1[,2:4]

outCC2pop <- read.csv("./boot data/outCC2pop_2018-05-15.csv")
outCC2pop <- outCC2pop[,2:4]
outCC2poplr <- read.csv("./boot data/outCC2poplr_2018-05-15.csv")
outCC2poplr <- outCC2poplr[,2:4]
outEG2pop <- read.csv("./boot data/outEG2pop_2018-05-15.csv")
outEG2pop <- outEG2pop[,2]
outEG2poplr <- read.csv("./boot data/outEG2poplr_2018-05-15.csv")
outEG2poplr <- outEG2poplr[,2]
outEP2pop <- read.csv("./boot data/outEP2pop_2018-05-15.csv")
outEP2pop <- outEP2pop[,2:4]
outEP2poplr <- read.csv("./boot data/outEP2poplr_2018-05-15.csv")
outEP2poplr <- outEP2poplr[,2:4]

sum_boot <- quantile(outCC1[,1], c(0.05,.5, 0.95))
sum_boot <- rbind(sum_boot, quantile(outCC1[,2], c(0.05,.5, 0.95)))
sum_boot <- rbind(sum_boot, quantile(outCC1[,3], c(0.05,.5, 0.95)))

sum_boot <- rbind(sum_boot, quantile(outEP1[,1], c(0.05,.5, 0.95)))
sum_boot <- rbind(sum_boot, quantile(outEP1[,2], c(0.05,.5, 0.95)))
sum_boot <- rbind(sum_boot, quantile(outEP1[,3], c(0.05,.5, 0.95)))

sum_boot <- rbind(sum_boot, quantile(outEG1, c(0.05,.5, 0.95)))

sum_boot <- rbind(sum_boot, quantile(outCC2pop[,1], c(0.05,.5, 0.95)))
sum_boot <- rbind(sum_boot, quantile(outCC2pop[,2], c(0.05,.5, 0.95)))
sum_boot <- rbind(sum_boot, quantile(outCC2pop[,3], c(0.05,.5, 0.95)))

sum_boot <- rbind(sum_boot, quantile(outEP2pop[,1], c(0.05,.5, 0.95)))
sum_boot <- rbind(sum_boot, quantile(outEP2pop[,2], c(0.05,.5, 0.95)))
sum_boot <- rbind(sum_boot, quantile(outEP2pop[,3], c(0.05,.5, 0.95)))

sum_boot <- rbind(sum_boot, quantile(outEG2pop, c(0.05,.5, 0.95)))

sum_boot <- rbind(sum_boot, quantile(outCC2poplr[,1], c(0.05,.5, 0.95)))
sum_boot <- rbind(sum_boot, quantile(outCC2poplr[,2], c(0.05,.5, 0.95)))
sum_boot <- rbind(sum_boot, quantile(outCC2poplr[,3], c(0.05,.5, 0.95)))

sum_boot <- rbind(sum_boot, quantile(outEP2poplr[,1], c(0.05,.5, 0.95)))
sum_boot <- rbind(sum_boot, quantile(outEP2poplr[,2], c(0.05,.5, 0.95)))
sum_boot <- rbind(sum_boot, quantile(outEP2poplr[,3], c(0.05,.5, 0.95)))

sum_boot <- rbind(sum_boot, quantile(outEG2poplr, c(0.05,.5, 0.95)))

sum_boot <- data.frame(matrix(sum_boot, ncol = 3))
sum_boot <- rename(sum_boot, c("X1"="5%", "X2"="Median", "X3"="95%"))
sum_boot$Model <- rep(c(rep("Climate Scepticsm", 3), rep("Environmental Protection", 3), rep("Environment vs. Growth", 1)),3)
sum_boot$'Independent Variable' <- c(rep("Populism", 14), rep("Populism x Political Ideology", 7))
sum_boot$'Dependent Variable' <- rep(c("not changing due to human activity", "is not changing", "Do not know", "Do not know", "Gone too far", "Not enough", "Trade-Off"), 3)
sum_boot$BES <- c(coef(mcc.1)[1,"populism"], coef(mcc.1)[2,"populism"], coef(mcc.1)[3,"populism"],
                  coef(meP.1)[1,"populism"], coef(meP.1)[2,"populism"], coef(meP.1)[3,"populism"],
                  coef(meG.1)["populism"],
                  coef(mcc.2)[1,"populism"], coef(mcc.2)[2,"populism"], coef(mcc.2)[3,"populism"],
                  coef(meP.3)[1,"populism"], coef(meP.3)[2,"populism"], coef(meP.3)[3,"populism"],
                  coef(meG.2)["populism"],
                  coef(mcc.2)[1,"populism:lr"], coef(mcc.2)[2,"populism:lr"], coef(mcc.2)[3,"populism:lr"],
                  coef(meP.3)[1,"populism:lr"], coef(meP.3)[2,"populism:lr"], coef(meP.3)[3,"populism:lr"],
                  coef(meG.2)["populism:lr"])

sum_boot <- data.frame(sum_boot$Model, sum_boot$'Independent Variable', sum_boot$'Dependent Variable', sum_boot$`5%`, sum_boot$Median, sum_boot$`95%`, sum_boot$BES)
sum_boot <- rename(sum_boot, c("sum_boot.Model"="Model", "sum_boot..Independent.Variable."="Independent Variable", "sum_boot..Dependent.Variable."="Dependent Variable", "sum_boot..5.."="5%", "sum_boot.Median"="Median", "sum_boot..95.."="95%", "sum_boot.BES" = "BES Coef."))

x.side <- xtable(sum_boot, caption = "Summary Bootstrapping", label = "sumBoot")
print(x.side,type = "html", file = "C:/Users/robhuber/polybox/Research/Populism Attitudes/Populism and Climate Scepticism/Submissions/PSRM/tables/boot.html", floating = TRUE, floating.environment = "sidewaystable")
print(x.side, floating = TRUE, floating.environment = "sidewaystable")

#Output ####

# CC1 ####

c <- c(outCC1$X1, outCC1$X2, outCC1$X3)
dv <- rep(c("is changing but not due to human activity", "is not changing", "Do not know"), each=nboot)
dv <- factor(dv, levels = c("is changing but not due to human activity", "is not changing", "Do not know"))

df_bootCC1 <- data.frame(c,dv)

c <- coef(mcc.1)[,"populism"]
dv <- factor(levels(dv))
txt <- "BES Coefficient"
df_hline <- data.frame(dv, c, txt)


p_bootCC1 <- ggplot(df_bootCC1, aes(x=c)) +
  geom_histogram(bins = 100)+
  theme_tufte() +
  geom_vline(xintercept=0, lty="dashed") +
  geom_vline(data = df_hline, aes(xintercept = c), color = "red", lwd=2) +
  ggtitle("On the subject of climate change do you think the world's climate ...") +
  theme(strip.background = element_rect(fill = 'white')) +
  labs(y="Count", x= "Coefficients") +
#  geom_text(data = df_hline, aes(label = txt, y= 60, x= c+0.025), color = "red", fontface = "bold", angle = 90) + 
  facet_grid(.~dv) +
  theme(text = element_text(size=20))
p_bootCC1

#CC2####

c <- c(outCC2pop[,1], outCC2pop[,2], outCC2pop[,3], outCC2poplr[,1], outCC2poplr[,2], outCC2poplr[,3])
dv <- rep(rep(c("is changing but not due to human activity", "is not changing", "Do not know"), each=nboot),2)
dv <- factor(dv, levels = c("is changing but not due to human activity", "is not changing", "Do not know"))
iv <- rep(rep(c("Populism", "Populism x Political Ideology"), each=3*nboot))
iv <- factor(iv, levels = c("Populism", "Populism x Political Ideology"))
df_bootCC2 <- data.frame(c,dv, iv)

c <- c(coef(mcc.2)[,c("populism")], coef(mcc.2)[,c("populism:lr")])
dv <- rep(factor(levels(dv), levels = c("is changing but not due to human activity", "is not changing", "Do not know")), 2)
iv <- rep(factor(levels(iv), levels = c("Populism", "Populism x Political Ideology")), each = 3)
txt <- "BES Coefficient"
df_hline <- data.frame(dv, iv, c, txt)


p_bootCC2 <- ggplot(df_bootCC2, aes(x=c))+
  geom_histogram(bins = 100)+
  facet_grid(dv~iv, scales = "free") +
  theme_tufte() +
  geom_vline(xintercept=0, lty="dashed") +
  geom_vline(data = df_hline, aes(xintercept = c), color = "red", lwd=2) +
  ggtitle("On the subject of climate change do you think the world's climate ...") +
  theme(strip.background = element_rect(fill = 'white')) +
  labs(y="Count", x= "Coefficients")+
  theme(text = element_text(size=20))
#  geom_text(data = df_hline, aes(label = txt, y=15, x= c+0.025), color = "red", fontface = "bold", angle = 90) 
p_bootCC2

# EP1 ####

c <- c(outEP1[,1], outEP1[,2], outEP1[,3])
dv <- rep(c("Do not know", "Gone too far", "Not enough"), each=nboot)
dv <- factor(dv, levels = c("Do not know", "Gone too far", "Not enough"))

df_bootEP1 <- data.frame(c,dv)

c <- coef(meP.1)[,"populism"]
dv <- factor(levels(dv))
txt <- "BES Coefficient"
df_hline <- data.frame(dv, c, txt)


p_bootEP1 <- ggplot(df_bootEP1, aes(x=c))+
  geom_histogram(bins = 100)+
  theme_tufte() +
  geom_vline(xintercept=0, lty="dashed") +
  geom_vline(data = df_hline, aes(xintercept = c), color = "red", lwd=2) +
  ggtitle("Do you think that measures to protect the environment have gone too far or not far enough?") +
  theme(strip.background = element_rect(fill = 'white')) +
  labs(y="Count", x= "Coefficients") +
#  geom_text(data = df_hline, aes(label = txt, y= 60, x= c+0.025), color = "red", fontface = "bold", angle = 90) + 
  #  geom_point(aes(y=60, x=mean(c)), size = 2) +
  #  geom_errorbarh(aes(xmin = as.numeric(quantile(c, c(0.025))), xmax = as.numeric(quantile(c, c(0.975))), height = 0, y=60), lwd = 2) +
  facet_grid(.~dv)+
  theme(text = element_text(size=20))
p_bootEP1

#EP2####

c <- c(outEP2pop[,1], outEP2pop[,2], outEP2pop[,3], outEP2poplr[,1], outEP2poplr[,2], outEP2poplr[,3])
dv <- rep(rep(c("Do not know", "Gone too far", "Not enough"), each=nboot), 2)
dv <- factor(dv, levels = c("Do not know", "Gone too far", "Not enough"))
iv <- rep(c("Populism", "Populism x Political Ideology"), each=3*nboot)
iv <- factor(iv, levels = c("Populism", "Populism x Political Ideology"))
df_bootEP2 <- data.frame(c,dv, iv)

#rm(df_hline)
c <- c(coef(meP.3)[,c("populism")], coef(meP.3)[,c("populism:lr")])
dv <- rep(factor(levels(dv), levels = c("Do not know", "Gone too far", "Not enough")), 2)
iv <- rep(factor(levels(iv), levels = c("Populism", "Populism x Political Ideology")), each = 3)
txt <- "BES Coefficient"
df_hline <- data.frame(dv, iv, c, txt)

p_bootEP2 <- ggplot(df_bootEP2, aes(x=c))+
  geom_histogram(bins = 100)+
  facet_grid(dv~iv, scales = "free") +
  theme_tufte() +
  geom_vline(xintercept=0, lty="dashed") +
  geom_vline(data = df_hline, aes(xintercept = c), color = "red", lwd=2) +
  ggtitle("Do you think that measures to protect the environment have gone too far or not far enough?") +
  theme(strip.background = element_rect(fill = 'white')) +
  labs(y="Count", x= "Coefficients")+
  theme(text = element_text(size=20))
#  geom_text(data = df_hline, aes(label = txt, y= 60, x= c+0.025), color = "red", fontface = "bold", angle = 90) 
p_bootEP2

# EG1 ####

c <- c(outEG1)
#dv <- rep(c("Do not know", "Gone too far", "Not enough"), each=nboot)
#dv <- factor(dv, levels = c("Do not know", "Gone too far", "Not enough"))

df_bootEG1 <- data.frame(c)

c <- coef(meG.1)["populism"]
#dv <- factor(levels(dv))
txt <- "BES Coefficient"
df_hline <- data.frame(c, txt)


p_bootEG1 <- ggplot(df_bootEG1, aes(x=c))+
  geom_histogram(bins = 100)+
  theme_tufte() +
  geom_vline(xintercept=0, lty="dashed") +
  geom_vline(data = df_hline, aes(xintercept = c), color = "red", lwd=2) +
  ggtitle("Environmental protection (0) vs. economic growth (10)") +
  theme(strip.background = element_rect(fill = 'white')) +
  labs(y="Count", x= "Coefficients")+
  theme(text = element_text(size=20))
#  geom_text(data = df_hline, aes(label = txt, y= 30, x= c+0.025), color = "red", fontface = "bold", angle = 90) 
p_bootEG1

#EG2####

c <- c(outEG2pop, outEG2poplr)
#dv <- rep(rep(c("Do not know", "Gone too far", "Not enough"), each=nboot), 2)
#dv <- factor(dv, levels = c("Do not know", "Gone too far", "Not enough"))
iv <- rep(c("Populism", "Populism x Political Ideology"), each=nboot)
iv <- factor(iv, levels = c("Populism", "Populism x Political Ideology"))
df_bootEG2 <- data.frame(c, iv)

#rm(df_hline)
c <- c(coef(meG.2)[("populism")], coef(meG.2)[("populism:lr")])
#dv <- rep(factor(levels(dv), levels = c("Do not know", "Gone too far", "Not enough")), 2)
iv <- rep(factor(levels(iv), levels = c("Populism", "Populism x Political Ideology")))
txt <- "BES Coefficient"
df_hline <- data.frame(iv, c, txt)

p_bootEG2 <- ggplot(df_bootEG2, aes(x=c))+
  geom_histogram(bins = 100)+
  facet_grid(~iv, scales="free") +
  theme_tufte() +
  geom_vline(xintercept=0, lty="dashed") +
  geom_vline(data = df_hline, aes(xintercept = c), color = "red", lwd=2) +
  ggtitle("Environmental protection (0) vs. economic growth (10)") +
  theme(strip.background = element_rect(fill = 'white')) +
  labs(y="Count", x= "Coefficients")+
  theme(text = element_text(size=20))
#  geom_text(data = df_hline, aes(label = txt, y= 30, x= c+0.035), color = "red", fontface = "bold", angle = 90) 
p_bootEG2

